One of the most annoying things one can encounter is being in the cinema, trying to enjoy a movie and having people around that talk. In this study we will try to link talking to the cinema with personality traits. Specifically, we asked participants to report how often they talk to the cinema in a 7 point likert scale and we administered questionnaires regarding agreeableness and consciousness which are part of the big-five perosnality traits. In addition, we administered questionnaires regarding the dark triad which refer to narcissism - Machiavellianism and psychopathy. We mesured also attention impulsivity which is related to psychopathy.

Our main hypotheses are:

Hypothesis 1: Narcissism will be related to talking to the cinema. Hypothesis 2: Psychopathy will be related to talking to cinema via a attention impulsivity Hypothesis 3: Agreeableness and consciousness will be relate to talking to the cinema.

Load libraries, set paths, load data

Libraries and paths


# Libraries
library(psych)      # for descriptive statistics and exploratory factor analysis (also needs *GPAroation* installed, will load it automatically, but not give error if it's not installed; so make sure you install this together with dependencies)
library(GPArotation)
library(dplyr)      # for data handling
library(data.table) # for re-arranging outputs of factor analysis to make nice plots (could be done in dplyr now, but not easily when originally coded)
library(ggplot2)    # for plotting
library(ggpubr)     # for plotting histograms
library(PerformanceAnalytics) # for plotting correlations
library(ggcorrplot) # for plotting correlations
library(outliers)   # different test to detect outliers
library(lavaan)     # for confirmatory factor analysis
library(ltm)        # for cronbach a
library(performance) # for inter-item correlation
library(lavaanPlot) # for plotting lavaan objects
library(semTools)   # for measurement invariance
library(brms)       # for bayesian regression



# Set paths
plotFolder                     = file.path(basepath,'/results')
statsFolder                    = file.path(basepath,'/figures')

Read file

#read csv file with raw data
setwd("~/Documents/cinema")
df <- read.csv(file="cinemaDat.csv", header=TRUE, sep=";",dec=",")

#Mahalobi’s distance - detecting outliers First check for multivariate outliers using mahalanobi’s distance.
The Mahalanobis distance is the distance between two points in a multivariate space. It’s often used to find outliers in statistical analyses that involve several variables. To determine if any of the distances are statistically significant, we need to calculate their p-values. The p-value for each distance is calculated as the p-value that corresponds to the Chi-Square statistic of the Mahalanobis distance with k-1 degrees of freedom, where k = number of variables. Typically a p-value that is less than .001 is considered to be an outlier.

mahalanobis distance

my.data<-df %>% dplyr::select(ATTENTION, EXTR, AGR, CON, MAC, NARC, PSYCH, AffEmp, CoEmp)

#check cor with outliers


#create new column in data frame to hold Mahalanobis distances
my.data$mahal <- mahalanobis(my.data, colMeans(my.data), cov(my.data))
#create new column in data frame to hold p-value for each Mahalanobis distance 
# df= k-1, k= number of variables that you use for testing 
my.data$p <- pchisq(my.data$mahal, df=8, lower.tail=FALSE) 

#for 8 df the critical value is 26.13

#drop outliers
df<-slice(df, -c(299, 257, 329,10, 131, 47))

In our sample, six participants can be considered as multivariate outliers according to mahalabis’ distance!

Plots and Correlations

Normal Distribution

Next we will proceed by checking if our data follow the normal distribution. We will do this by checking skeweness and kirtosis and by also checking the qq plots for each variables. For sample sizes greater than 300, depend on the absolute values of skewness and kurtosis without considering z-values. Either an absolute skew value larger than 2 or an absolute kurtosis (proper) larger than 7 may be used as reference values for determining substantial non-normality. For more information see: Statistical notes for clinical researchers: assessing normal distribution (2) using skewness and kurtosis (Kim, 2013)

Correlations

Next we will check the correlations. As our depended variable (talking to cinema) is ordinal we will use sperman sho as method. We will also plot the correlation, scatter plots and histograms of the variables of intererest.

Descriptives statistics, plots and correlations


my.data<-df %>% dplyr::select(food, cell, Talk, ATTENTION, EXTR, AGR, CON, MAC, NARC, PSYCH, AffEmp, CoEmp)
describe(my.data)

#plot qq plots for normality
qq.plots<-list()

qq.plots$ATTENTION<-ggqqplot(df$ATTENTION, ylab = "Attentional Impulsiveness")
qq.plots$EXTR<-ggqqplot(df$EXTR, ylab = "Extraversion")
qq.plots$CON<-ggqqplot(df$CON, ylab = "consientousness")
qq.plots$AGR<-ggqqplot(df$AGR, ylab = "Agreeableness")
qq.plots$MAC<-ggqqplot(df$MAC, ylab = "Maciavelism")
qq.plots$NARC<-ggqqplot(df$NARC, ylab = "Narcissism")
qq.plots$PSYCH<-ggqqplot(df$PSYCH, ylab = "Psychopathy")
qq.plots$AffEmp<-ggqqplot(df$AffEmp, ylab = "Affective Empathy")
qq.plots$CoEmp<-ggqqplot(df$CoEmp, ylab = "Cognitive Empathy")

plot.normality<-ggarrange(qq.plots$ATTENTION,qq.plots$EXTR,qq.plots$CON,qq.plots$AGR,qq.plots$MAC,qq.plots$NARC, qq.plots$PSYCH,qq.plots$AffEmp, qq.plots$CoEmp, nrow = 3,ncol=3, hjust = 0)
plot(plot.normality)



#cordata
cordata<-(my.data)

#plot correlation
chart.Correlation(cordata, method =  "spearman", histogram = TRUE, pch=19)


cordata<-(my.data)

Hypotheses Testing

In the next step we will check our hypotheses.

Hypothesis 1.

The first hypothesis states that narcissism will be related to talking in cinema. To test this hypothesis we will use Bayesian regression. We will use talking as depended variable and narcissism as the regressor. Moreover, we will control for covoriates of no interest such as age, gender and education. We will use education as monotonic predictor. Since talking is an ordinal variable we will use family = “cumulative”. We will assess if the one-tailed 95% Bayesian credible intervals of the regression weights exclude zero (significant) or include zero (not significant).



# Settings for BRMS
niter=8000
nchains=4
adaptdeltas=0.9




library(brms)
h1.<-brm(formula=Talk
                   ~1 + NARC+gender+age+mo(educ), data=df,save_all_pars = TRUE, family='cumulative',iter=niter,chains=nchains,cores=nchains,control=list(adapt_delta=adaptdeltas))
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include   -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1

summary(h1.)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Talk ~ 1 + NARC + gender + age + mo(educ) 
##    Data: df (Number of observations: 352) 
##   Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
##          total post-warmup draws = 16000
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]    -0.35      0.82    -1.90     1.30 1.00    10261     9701
## Intercept[2]     1.18      0.82    -0.38     2.83 1.00    10260    10070
## Intercept[3]     1.95      0.82     0.38     3.61 1.00    10251     9510
## Intercept[4]     2.86      0.83     1.28     4.56 1.00    10138     9788
## Intercept[5]     4.01      0.86     2.38     5.75 1.00     9831     8948
## Intercept[6]     5.49      0.95     3.70     7.46 1.00     9788     9120
## NARC             0.35      0.11     0.13     0.57 1.00    14181    10504
## gender           0.27      0.27    -0.26     0.81 1.00    12636    10560
## age             -0.03      0.01    -0.05    -0.01 1.00    12310    11912
## moeduc           0.43      0.44    -0.16     1.41 1.00     2395     5430
## 
## Simplex Parameters: 
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## moeduc1[1]     0.13      0.13     0.00     0.47 1.00     4425     3085
## moeduc1[2]     0.07      0.10     0.00     0.39 1.00     3929     8358
## moeduc1[3]     0.07      0.10     0.00     0.37 1.00     4220     7892
## moeduc1[4]     0.08      0.10     0.00     0.37 1.00     5373     8389
## moeduc1[5]     0.37      0.23     0.01     0.80 1.00     3336     3031
## moeduc1[6]     0.27      0.20     0.01     0.72 1.00     7221     8130
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Hypothesis 2.

Hypothesis 2 states that psychopathy is related to talking via attentions impulsivity. To test this hypothesis we will use psychopathy as a predictor, talking as independent variable and attention implulsivity as mediator. We will treat the afforementioned variables as latent variables extracted by the observed variables (the items of each questionnaire). Since both the items as well as the dependend variable are ordinal we will use DWLS estimator which is suitable for ordinal data, (see Mîndrilă, 2010; Cheng-Hsien Li, 2016)

To assess good-fitting of the model χ2 and its degrees of freedom (df) were used. For χ2 values associated with p>0,5 were considered good-fitting models, although it has to be mentioned that the p value of this test is sensitive to large sample size. In addition, the root mean square error of approximation (RMSEA) with its 90% confidence intervals (CI), the standardized root mean square residuals (SRMR) and the comparative fit index (CFI) were used. For RMSEA and SRMR values values < 0.08 are acceptable. For CFI values > 0.90 were considered as indicators of good fit (Brown, 2006).

Regarding the mediation, we computed the cross product of the two direct paths coefficient to obtain the indirect path coefficient. Statistical significance level was set at α=0.05.

model2<-'
         psych=~Psych1+Psych2_r+Psych3+Psych4+Psych5+Psych6+Psych7_r+Psych8+Psych9
         


Attentio=~Att1+Att2+Att3_r+Att4+Att5_r+Att6+Att7+Att8
       Attentio~a*psych
    
         Talk~b*Attentio
         Talk~c*psych
       
         
          ab :=a*b
total :=c+(a*b)

'

sem.fit2 <- sem(model2, data=df ,estimator = "DWLS")
summary(sem.fit2, fit.measures=TRUE, standardized = TRUE,rsq=TRUE)
## lavaan 0.6-11 ended normally after 68 iterations
## 
##   Estimator                                       DWLS
##   Optimization method                           NLMINB
##   Number of model parameters                        38
##                                                       
##   Number of observations                           352
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                               200.146
##   Degrees of freedom                               133
##   P-value (Chi-square)                           0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                              1252.229
##   Degrees of freedom                               153
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.939
##   Tucker-Lewis Index (TLI)                       0.930
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.038
##   90 Percent confidence interval - lower         0.027
##   90 Percent confidence interval - upper         0.048
##   P-value RMSEA <= 0.05                          0.973
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.062
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   psych =~                                                              
##     Psych1            1.000                               0.854    0.499
##     Psych2_r         -0.456    0.085   -5.347    0.000   -0.390   -0.226
##     Psych3            0.994    0.120    8.268    0.000    0.849    0.625
##     Psych4            0.856    0.104    8.233    0.000    0.731    0.597
##     Psych5            1.490    0.173    8.604    0.000    1.273    0.666
##     Psych6            1.059    0.129    8.236    0.000    0.905    0.539
##     Psych7_r         -0.423    0.086   -4.901    0.000   -0.361   -0.199
##     Psych8            0.567    0.097    5.861    0.000    0.484    0.267
##     Psych9            1.115    0.134    8.348    0.000    0.953    0.650
##   Attentio =~                                                           
##     Att1              1.000                               0.215    0.253
##     Att2              1.204    0.263    4.574    0.000    0.259    0.269
##     Att3_r           -1.691    0.328   -5.157    0.000   -0.364   -0.452
##     Att4              2.692    0.499    5.389    0.000    0.579    0.618
##     Att5_r           -1.100    0.242   -4.544    0.000   -0.237   -0.300
##     Att6              1.238    0.273    4.533    0.000    0.266    0.302
##     Att7              1.931    0.372    5.191    0.000    0.415    0.444
##     Att8              3.082    0.567    5.433    0.000    0.663    0.781
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   Attentio ~                                                            
##     psych      (a)    0.080    0.017    4.779    0.000    0.319    0.319
##   Talk ~                                                                
##     Attentio   (b)    2.072    0.501    4.135    0.000    0.446    0.299
##     psych      (c)    0.018    0.078    0.226    0.821    0.015    0.010
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .Psych1            2.203    0.260    8.487    0.000    2.203    0.751
##    .Psych2_r          2.814    0.198   14.235    0.000    2.814    0.949
##    .Psych3            1.126    0.257    4.379    0.000    1.126    0.609
##    .Psych4            0.968    0.188    5.148    0.000    0.968    0.644
##    .Psych5            2.034    0.308    6.614    0.000    2.034    0.556
##    .Psych6            2.003    0.254    7.879    0.000    2.003    0.710
##    .Psych7_r          3.176    0.339    9.376    0.000    3.176    0.961
##    .Psych8            3.059    0.266   11.510    0.000    3.059    0.929
##    .Psych9            1.240    0.254    4.890    0.000    1.240    0.577
##    .Att1              0.675    0.058   11.601    0.000    0.675    0.936
##    .Att2              0.859    0.055   15.665    0.000    0.859    0.928
##    .Att3_r            0.517    0.050   10.242    0.000    0.517    0.796
##    .Att4              0.543    0.079    6.846    0.000    0.543    0.618
##    .Att5_r            0.565    0.043   13.078    0.000    0.565    0.910
##    .Att6              0.707    0.058   12.132    0.000    0.707    0.909
##    .Att7              0.702    0.061   11.440    0.000    0.702    0.803
##    .Att8              0.281    0.082    3.430    0.001    0.281    0.390
##    .Talk              2.013    0.169   11.918    0.000    2.013    0.908
##     psych             0.730    0.127    5.736    0.000    1.000    1.000
##    .Attentio          0.042    0.014    3.025    0.002    0.898    0.898
## 
## R-Square:
##                    Estimate
##     Psych1            0.249
##     Psych2_r          0.051
##     Psych3            0.391
##     Psych4            0.356
##     Psych5            0.444
##     Psych6            0.290
##     Psych7_r          0.039
##     Psych8            0.071
##     Psych9            0.423
##     Att1              0.064
##     Att2              0.072
##     Att3_r            0.204
##     Att4              0.382
##     Att5_r            0.090
##     Att6              0.091
##     Att7              0.197
##     Att8              0.610
##     Talk              0.092
##     Attentio          0.102
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     ab                0.167    0.038    4.418    0.000    0.142    0.096
##     total             0.184    0.065    2.815    0.005    0.157    0.106
lavaanPlot(model=sem.fit2, coefs = TRUE, stars = "regress", digits = 2, stand = TRUE)

Hypothesis 3

Hypothesis 3 states that Agreeableness and consciousness will be relate to talking to the cinema. To test this hypothesis we will do model comparison. Models give means of putting scientific hypotheses into practice; as a result, model comparison and hypothesis testing are intertwined (Bruno Nicenboim, Daniel Schad, and Shravan Vasishth, 2022).
The Bayes factor reflects a ratio that provides information about how much more likely the observed data are between two compared models. Bayes factor greater than 1 provides evidence in support of one model over the other. We will use Bayes factor to compare the model described in the above-mentioned hypothesis while controlling for covariates of no interesting, against a null model that involves only the intercept without any predictor and only covariates of no interested (age, gender and education treated as monotonic). For more information see vignettes/bayes_factors.Rmd https://easystats.github.io/bayestestR/articles/bayes_factors.html

Also for more information regarding model comparison see:

  1. Makowski D, Ben-Shachar MS, Lüdecke D. bayestestR: Describing Effects and their Uncertainty, Existence and Significance within the Bayesian Framework. J Open Source Softw. 2019 Aug 13;4(40):1541.
  2. Vehtari A, Gelman A, Gabry J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat Comput. 2017 Sep;27(5):1413–32.


#fit the null model
null.h3<-brm(formula=Talk
                   ~1 + +gender+age+mo(educ), data=df,save_all_pars = TRUE, family='cumulative',iter=niter,chains=nchains,cores=nchains,control=list(adapt_delta=adaptdeltas))
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include   -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
summary(null.h3)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Talk ~ 1 + +gender + age + mo(educ) 
##    Data: df (Number of observations: 352) 
##   Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
##          total post-warmup draws = 16000
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]    -1.57      0.71    -2.90    -0.12 1.00     6682     4240
## Intercept[2]    -0.08      0.70    -1.40     1.38 1.00     6566     4421
## Intercept[3]     0.68      0.70    -0.66     2.13 1.00     6657     4574
## Intercept[4]     1.58      0.71     0.23     3.05 1.00     6214     4403
## Intercept[5]     2.72      0.74     1.33     4.24 1.00     6335     4153
## Intercept[6]     4.19      0.85     2.60     5.94 1.00     6469     4653
## gender           0.22      0.27    -0.31     0.75 1.00    14921    11865
## age             -0.03      0.01    -0.05    -0.01 1.00    14143    10820
## moeduc           0.38      0.43    -0.16     1.38 1.00     2621     3770
## 
## Simplex Parameters: 
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## moeduc1[1]     0.14      0.13     0.00     0.48 1.00     7141     7745
## moeduc1[2]     0.08      0.11     0.00     0.40 1.00     4194     8314
## moeduc1[3]     0.08      0.10     0.00     0.39 1.00     4490     7935
## moeduc1[4]     0.09      0.10     0.00     0.39 1.00     6240     7770
## moeduc1[5]     0.34      0.22     0.01     0.77 1.00     4797     7732
## moeduc1[6]     0.27      0.21     0.01     0.73 1.00     6570     7433
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
        

#fit the h3 model
h3.<-brm(formula=Talk
                   ~1 + AGR+CON+gender+age+mo(educ), data=df,save_all_pars = TRUE, family='cumulative',iter=niter,chains=nchains,cores=nchains,control=list(adapt_delta=adaptdeltas))

summary(h3.)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Talk ~ 1 + AGR + CON + gender + age + mo(educ) 
##    Data: df (Number of observations: 352) 
##   Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
##          total post-warmup draws = 16000
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]    -3.92      1.00    -5.88    -1.93 1.00     9213     5060
## Intercept[2]    -2.39      0.99    -4.32    -0.41 1.00     9706     5261
## Intercept[3]    -1.61      0.98    -3.53     0.34 1.00    10070     5293
## Intercept[4]    -0.69      0.98    -2.59     1.27 1.00     9621     5302
## Intercept[5]     0.46      1.00    -1.45     2.47 1.00     9326     5109
## Intercept[6]     1.94      1.09    -0.13     4.18 1.00     7751     3005
## AGR             -0.21      0.12    -0.45     0.03 1.00    15085    10327
## CON             -0.30      0.11    -0.51    -0.09 1.00    12678    10741
## gender           0.20      0.27    -0.33     0.71 1.00    14251    10518
## age             -0.02      0.01    -0.04    -0.00 1.00    14470    10257
## moeduc           0.42      0.43    -0.10     1.43 1.00     3079     5146
## 
## Simplex Parameters: 
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## moeduc1[1]     0.14      0.13     0.00     0.47 1.00     6253     6638
## moeduc1[2]     0.07      0.10     0.00     0.36 1.00     5758     8898
## moeduc1[3]     0.07      0.09     0.00     0.34 1.00     6220     9164
## moeduc1[4]     0.09      0.10     0.00     0.38 1.00     7207     7549
## moeduc1[5]     0.34      0.21     0.02     0.76 1.00     6932     8011
## moeduc1[6]     0.29      0.21     0.01     0.75 1.00     6187     4948
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).



# Compute Bayes factor
bayes_factor(h3., null.h3)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6
## Iteration: 7
## Estimated Bayes factor in favor of h3. over null.h3: 38.96313